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2  Technical  accomplishments 

Granular  energetic  materials  exhibit  com¬ 
plex  chaotic  behavior  due  to  the  coexistence 
of  a  wide  range  of  energy  scales  without 
scale  separation.  The  main  challenges  in¬ 
volved  in  modeling  the  physical  processes 
leading  to  initiation  of  explosive  reactions 
are  (i)  the  lack  of  a  general  model  for  het¬ 
erogeneous  granular  media  under  compac¬ 
tion  and  (ii)  the  lack  of  a  reliable  multi¬ 
scale  discrete-to-continuum  framework  for 
describing  diffusion-advection-reaction  pro¬ 
cesses  in  complex  particulate  media.  Most 
conventional  methods  for  studying  visco¬ 
plastic  deformations  of  granular  media  un- 
,  ,  ,  .....  , .  Figure  1 :  Interaction  between  granular  dissipation 

der  shear  and  compression  in  ideal  condi-  & 

,  .  ..  ....  .  .  ,  -  .  due  to  friction,  stress  waves  and  heat  diffusion, 

tions  overlook  the  effect  of  spatial  hetero¬ 
geneity  in  granular  structure.  This  heterogeneity  is  believed  to  play  a  major  role  in  stress  and 
heat  localization  events  responsible  for  initiating  reactions  in  energetic  materials.  In  particu¬ 
lar,  it  has  been  observed  that  so-called  “hot-spots”  emerge  as  a  consequence  of  visco-plastic 
pore  collapse,  inter-granular  friction,  and  granular  compaction.  The  wide  range  in  stress, 
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strain,  and  dissipation  found  in  energetic  materials  magnify  the  multiscale  behavior  of  the 
system.  That  is,  microscopic  events  invariably  affect  the  macroscopic  behavior  of  the  sys¬ 
tem  and  cannot  be  neglected,  yet  are  also  impossible  to  predict  deterministically.  Specifically, 
macroscopic  stress  boundary  conditions  induce  heterogeneous  deformations  at  the  grain  level 
which  causes  friction  and  grain  deformations  at  the  microscopic  level.  This  generates  thermal 
fluctuations  and  chemical  reactions  at  the  molecular  level,  which  in  turn  builds  up  a  shock 
wave  whose  power  is  several  orders  of  magnitude  higher  than  the  initial  conditions. 

During  the  two  years  of  this  project  we  investigated  two  of  the  processes  in  Figure  1: 
multiscale  diffusion  /  heat  transfer  and  multiscale  dynamics  of  granular  materials. 


2.1  Hybrid  Discrete-Continuum  Models  of  Heat  Dissipation 

Efficient  coupling  of  continuum  (deterministic  or  stochastic)  constitutive  solvers  with  their 
discrete  (stochastic,  particle-based)  counterparts  is  a  common  challenge  in  hybrid  and  mul¬ 
tiphysics  simulations.  We  [Bakarji  and  Tartakovsky,  2016]  studied  interfacial,  tightly  cou¬ 
pled  simulations  of  diffusion  that  combine  continuum  and  particle-based  solvers.  The  latter 
employed  the  reverse  Brownian  motion  (rBm),  a  Monte  Carlo  approach  that  allows  one  to 
enforce  inhomogeneous  Dirichlet,  Neumann,  or  Robin  boundary  conditions  and  is  trivially 
parallelizable.  In  Brownian  motion,  a  particle’s  trajectory  X(t)  evolves  in  time  according  to 
a  stochastic  differential  equation  dX(i)  =  y/2 ad  dW(£)  where  ad  is  a  diffusion  coefficient, 
and  dW(t)  ~  A/*(0.  d t)  is  a  d-dimensional  Wiener  process.  Our  Monte  Carlo  simulations 
used  the  rBm  implementation,  in  which  individual  trajectories  of  NMc  particles  released  at 
point  x  at  time  t  satisfy  X(f  —  Atd)  =  X(f)  —  v/2c^A/’(0,  Atd).  Given  an  initial  condition 
U[n (x) ,  and  the  functions  uD(x,  t)  and  JN(x,  t)  prescribed  respectively  on  the  Dirichlet  and 
Neumann  boundary  conditions,  the  sample  mean  temperature  fid(x,  t)  at  space-time  point 
(x,  t)  is  computed  as  a  weighted  sum 


fid(x,t) 


N-m 


Ni 


'  + 


MC 


N&  c  iVN 

jv  +  ~rj—  On 

iVMC  JVmC 
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of  sample  averages  of  the  initial  and  boundary  functions,  uin[Xj(0)],  uD[X.l(t  —  Tt),  t  —  T] 
and  JN[Xj(f  —  Titj),t  —  Ti:j ].  Here  Nin,  ND  and  Nn  are  the  numbers  of  particles  that  reached 
the  initial  state  and  the  boundaries,  respectively;  and  7]  is  the  /th  particle’s  exit  time. 

We  developed  a  number  of  numerical  approaches  for  improving  the  accuracy  of  rBm 
in  the  presence  of  inhomogeneous  Neumann  boundary  condition  and  alternative  strategies 
for  coupling  the  rBm  solver  with  its  continuum  counterpart.  Numerical  experiments  were 
used  to  investigate  the  convergence,  stability,  and  computational  efficiency  of  the  proposed 
hybrid  algorithm.  Our  analysis  revealed  that  the  use  of  Monte  Carlo  simulations  based  on  the 
reverse  Brownian  motion  (rBm)  in  the  context  of  discrete-to-continuum  hybrid  simulations 
has  a  number  of  advantages.  These  include 


1.  the  ability  to  use  very  large  hybrid  time  step  A 4  without  compromising  accuracy, 

2.  the  ability  to  compute  the  solution  only  near  the  boundaries  to  ensure  the  continuity  of 
the  flux  in  the  hybrid  method, 
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3.  the  ability  to  use  the  continuum  domain  as  a  deterministic  source  of  Dirichlet,  Neumann 
and  initial  boundary  conditions  for  the  rBm,  and  therefore 

4.  a  controllable  loss  of  accuracy  given  a  flexible  choice  of  NMC  at  every  location  in  the 
particle  domain. 

Our  hybrid  algorithm  is  easy  to  implement  in  any  number  of  dimensions.  Furthermore,  ex¬ 
tending  the  hybrid  model  to  advection-diffusion  equations  is  relatively  straightforward. 

2.2  Hybrid  Discrete- Continuum  Models  of  Compaction  of  Granular  Ma¬ 
terials 

Given  a  force  applied  to  the  top  of  a  cylinder  filled  with  monodisperse  granular  medium 
with  spherical  metal  beads  (Figure  2),  find  the  height  of  the  top  lid  h(t)  as  a  function  of 
time.  The  particles  are  assumed  to  be  formed  of  a  deformable  and  incompressible  metal 
(e.g.,  aluminum).  The  pressure  P  on  the  top  is  assumed  to  be  constant.  It  is  expected  that 
the  presence  of  a  large  pore  will  induce  more  compactions,  i.e.,  hV0K  <  /ihom.  We  estimated 
the  height  of  the  lid  in  the  presence  of  a  macro-pore  as  compared  to  the  case  without  macro¬ 
pores.  More  specifically,  we  studied  the  effects  of  heterogeneity  by  using  a  local  model  for 
pore  collapse  in  a  granular  medium. 

Applied  Pressure 

III  n  III  III  n  111 


Skorokhod- 

Olevsky 

Model 


Macropore 


Figure  2:  A  physical  system  (left)  and  its  hybrid  model  (right). 

A  discrete  representation  of  a  granular  material’s  dynamics  was  based  on  the  Carroll-Holt 
(CH)  model,  which  is  an  axisymmetric  elasto-plastic  model  of  the  collapse  of  an  incompress¬ 
ible  metallic  shell.  We  used  the  modified  CH  model  that  introduces  a  restriction  on  the  extent 
to  which  the  pore  can  collapse,  avoiding  singularity  at  its  center.  While  using  this  model,  the 
coupling  has  to  be  done  in  a  way  that  transforms  the  compressible  shell  of  the  actual  case 
into  an  incompressible  shell. 
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Figure  3:  Coupling  algorithm  for  hybrid  simulations. 


We  used  the  non-linear  viscous  model  of  Skorokhod  and  Olevsky  as  a  continuum  (macro¬ 
scopic)  description  of  granular  compaction.  The  model  posits  that  stress  (al3)  depends  on 
strain  rate  (e^),  and  that  bulk  (7)  and  shear  (</>)  moduli  depend  on  porosity  (0): 


aij  ~ 


a(W) 

W~ 


(f>(9)eij  +  ( (f)(6)  —  -^ip(6)  )  eSy 


e  = 


duj 

dx. 


where  is  the  Kronecker  delta  function;  and 


4>  =  <j>8(  i  -  ef 


v  3  0 


(2) 


(3) 


In  the  linear- viscous  case,  a(W)  =  r/0W;  in  the  perfectly  plastic  case,  a(W)  =  ay  where 
W  =  y/(l  —  #)/[(/>(#)72  +  ■0(6,)e2]  and  7  =  '/¥!'.  The  continuity  equation  reduces  to  e  = 
0/(1  —  0).  These  constitutive  equations  were  then  used  in  the  Cauchy  equation  of  motion. 

These  two  levels  of  description  were  combined  in  a  globally  energy-conserving  hybrid 
model,  whose  schematic  representation  is  provided  in  Figure  3. 


Figure  4:  (a)  Changes  in  porosity  predicted  with  a  superposition  of  Carroll-Holt  solutions, 
(b)  Piston’s  position  computed  with  the  hybrid  model. 

This  hybrid  model  was  first  used  to  simulate  compaction  of  a  granular  material  with 
incompressible  matrix.  In  this  case,  compression  is  only  due  to  the  macropores  and  the 
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stress  distribution  around  the  pores  can  be  deduced  by  assuming  the  pores  to  be  in  an  infinite 
medium  (compared  to  the  whole  domain).  Figure  4a  shows  the  combined  effect  of  multiple 
macropores  using  a  superposition  of  Carroll-Holt  models  with  different  radii  and  surface 
stresses.  This  gives  an  idea  on  what  to  expect  in  the  difference  between  a  homogeneous  and 
heterogeneous  compaction. 

Assuming  a  linear- viscous  continuum  model  for  the  matrix,  we  explored  the  dependence 
of  the  height  of  the  lid  h(t)  as  a  function  of  viscosity.  Figure  4b  shows  a  range  of  viscosity 
of  one  order  of  magnitude. 


2.3  Fluctuating  Macroscopic  Models 

Inspired  by  fluctuating  Navier-Stokes  equations  of  hydrodynamics,  we  explored  ways  to 
model  unresolved  micro-scales,  e.g.,  heterogeneity  due  to  macropores,  as  a  random  source  in 
the  Cauchy  equation  of  motion, 

=  V-cr+2pf(x,f)  (4) 

where  p  is  the  density,  u  is  the  velocity  at  point  x  and  time  t,  and  the  stress  tensor  cr  is 
related  to  the  strain  rate  e  by  a  constitutive  law,  e.g.,  by  the  Skorokhod-Olevsky  relation  (2). 
The  (possibly  random)  indicator  function  2p(x)  for  a  macropore  region  f2p(t)  is  defined  as 
2p  =  1  if  x  £  f2p  and  =  0  otherwise.  The  macropore  region  Op(t)  is  a  multi-connected 
domain  comprising  multiple  macropores,  which  can  be  randomly  distributed  throughout  the 
material.  The  random  source  vector  f(x,  t)  is  treated  as  zero-mean  white  noise,  E[f  (x,  t)]  = 
0,  E[/j(x,  t)  fi( y,  t)]  =  vfS(x  —  y )S(t  —  r)  where  vf  is  the  variance  of  the  ith  component  of 
the  noise  and  5(  )  is  the  Dirac  delta  function.  Figure  5  exhibits  an  average  velocity  distribution 
within  the  granular  material  undergoing  slow  compaction,  for  the  case  of  a  single  macropore 
and  given  noise  strength  vf. 

111.111 


Figure  5:  Average  velocity  distribution  within  the  granular  material  undergoing  compaction. 

Our  still  elusive  goal  is  to  relate  the  variance  vf  to  material  properties,  e.g.,  grain-size 
distribution,  via  the  fluctuation-dissipation  theorem. 
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